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We study fractional Brownian motion (fBm) characterized by the Hurst exponent H. Using 
a Monte Carlo sampling technique, we are able to numerically generate fBm processes with an 
absorbing boundary at the origin at discrete times for a large number of 10 7 time steps even for 
small values like H = 1/4. The results are compatible with previous analytical results that the 
distribution of (rescaled) endpoints y follow a power law P+(y) ~ / with <p = (l-H)/H, even for 
small values of H. Furthermore, for the case H = 0.5 we also study analytically the finite-length 
corrections to the first order, namely a plateau of P+ (y) for y — > which decreases with increasing 
process length. These corrections are compatible with the numerical results. 
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I. INTRODUCTION 

The Brownian motion plays a key role in modern the- 
oretical physics, as it explains many effects observed in 
physical systems. It is currently used in various fields of 
science to understand, for instance, the trend of financial 
markets, the dynamics of complex molecules within the 
cell or the animals food-searching strategies. In order to 
describe the fluctuations in these systems, it is often nec- 
essary to go beyond the Brownian motion and consider 
random walkers whose mean square displacement grows 
over time in a nonlinear way. The term used to refer 
to this situation is anomalous diffusion, in particular sub 
diffusion if the mean square displacement grows less than 
linearly, super diffusion if it is faster. 

In practice, anomalous diffusion occurs whenever the 
process x(t) is self affine (at least at large time) with 
a characteristic value of the so-called Hurst exponent 
H 7^ 1/2, so that the displacement grows with time as 
t h . A remarkable example of process displaying anoma- 
lous diffusion is the fractional Brownian motion (fBm), 
originally introduced by Mandelbrot This process is 
self-affine Gaussian process with < H < 1. A Gaus- 
sian process is completely defined by its autocorrelation 
function, which for fBm writes as 



(x(r)x(r')) 



2H 



+ (r') 



2H 



J\1B 



(1) 



where x(0) = and the brackets (...) refer to an en- 
semble average over the realizations of the Gaussian 
processes. The strength of the correlation is described 
by the Hurst exponent. Note that Eq. {T]) implies 
([ar(Ti) - x(t 2 )} 2 ) = 2|n - t 2 \ 2H . This means, H = 1/2 
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corresponds to the Brownian motion (standard diffu- 
sion), H > 1/2 to super-diffusive paths, while H < 1/2 
correspond to sub-diffusive paths. 

Recently these random walks have found to be relevant 
for many physical applications. Among them we mention 
the fluctuations of a tagged monomer of an equilibrated 
Rouse chain 0, 0] or of a tagged particle in the one di- 
mensional system [1, Q . In both cases the motion of the 
tagged object can be modeled as a fractional Brownian 
motion with H — 1/4. Other physical processes such as 
the mechanical unzipping of DNA , the translocation of 
polymers through nanopores @, [7H9| , and subdiffusion of 
macromolecules inside cells and membranes [l(3-[l2j can 
be well described by fBm diffusion. 

Despite the large number of cases where fBm is ob- 
served, very little is known about the properties of this 
process when is confined in a domain of the space, which 
is the case for many of the above mentioned applica- 
tions. A concrete example is given by the process of poly- 
mer translocation. In this case it has been shown [2J, [7| 
that the fraction of the polymer penetrated inside the 
nanoporc fluctuates with time as a fBm walker confined 
inside the interval (which corresponds to the translo- 
cation failure) and 1 (which corresponds to completed 
translocation). 

In presence of boundaries the translational invariance 
is lost and analytical representations like the fractional 
Langevin equation are of little help [l3[ . In these situa- 
tions, numerical simulations remain the option to answer 
many concrete questions arising from biology and physics 
[b4T - [l6| . Recently, methods where developed to stdy fBm 
for a system with a potential [TtJ and for fBm in con- 
fined geometries (i~8| . Here we present a new generic 
numerical method which we use here to study these pro- 
cesses in presence of boundaries. We will study in detail 
the case where there is an absorbing wall in x = 0: we 
thus consider only the paths that remain positive up to 
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the final time r. Recently, analytical predictions [7|, |19( 
were obtained for the distribution P+{y) of rescaled mo- 
tion endpoints y ~ x/(t h ) at end time r. A possible 
numerical strategy consists in the direct sampling of in- 
step fBm paths xq,x\, . . . ,xl at discrete times, starting 
at xo = 0. This strategy is demanding, in particular 
for H < 1/2, since in presence of an absorbing bound- 
ary the success probability of generating a non-absorbed 
trajectory is very small. Hence, such simulations were 
restricted to a small number L of discrete steps. Here, 
using a Markov Chain approach, we were able to gener- 
ate long fBm processes up to L = 10 7 discrete steps for 
values such as H = 1/4, 2/5, 1/2, 2/3 and 3/4. 

The outline of the remainer of the paper is as follows: 
Next, we present our numerical approach and then we 
present our numerical results: First - in order to verify 
that our approach is working- we consider the case of 
Brownian motion (H — 1/2), where analytical results for 
the finite-length corrections are available. Furthermore 
we study the cases H = 1/4, H = 2/5 and H = 2/3, 
H = 3/4 as examples for the two regions H < 1/2 and 
H > 1/2. In all cases the results are compatible with 
previous analyitcial predictions. Finally, we summarize 
our results. 



II. NUMERICAL METHODS 

To generate fBm processes on a computer, we study 
discrete-time random walks with suitable correlations. It 
is useful to introduce the increments of the random walk, 
namely Axi — xj+i — xi. For Gaussian processes the 
increments are Gaussian variables defined by their auto- 
correlation function. Using Eq. (TTJ we can compute the 
autocorrelation function of the Gaussian increments: 

Ci +m j = {Axi +m Ax{) (2) 
= |m + 1\ 2H - 2\m\ 2H + \m- 1\ 2H = C(m) , 

We note that this function is independent of the initial 
time I. Matrices having this property are called Toeplitz 
matrices. Moreover thus the increments are identical 
Gaussian numbers with variance a 2 = 2 displaying power 
law correlations. Taking the limit m — > oo it is easy to 
extract the power law decay of these correlations. For 
super diffusive fBm (H > 1/2), C(m) is positively cor- 
related with a decay as m^ 2 ^ H \ Positive correlation 
means that there is a high probability to observe a long 
sequence of increments of same sign. For sub diffusive 
fBm (H < 1/2), C(m) is negatively correlated and decay 
as —mr 2 ( x ~ HS >. Negative correlation means that there is 
a high probability to observe a long sequence of incre- 
ments of oscillating sign. 

The direct generation of L steps with increment corre- 
lation Eq.Q is straightforward, in principle. The start- 
ing point is a vector £ = (£o, £i, • • • , of L inde- 
pendent and identically distributed (iid) Gaussian (mean 
zero, variance one) numbers ~ G(0, 1). For the uncorre- 



cted case (H = 1/2) one could directly use the ran- 
dom numbers, multiplied by \[2 to obtain the right Ci : i, 
as increments of the fBm processes, i.e., x unco "(L) — 

For the case H ^ 1/2, since C as a correlation ma- 
trix is positive semi-definite, there exist a matrix A such 
that C = A 2 . Thus, one could use Ax = A£ to obtain a 
random vector with the desired property Eq. @. Nev- 
ertheless, this is too time consuming, since it requires 
diagonalizing a L x L C matrix once (~ L 3 operations) 
and, for each process, the multiplication with the L x L 
matrix, A (~ L 2 operations) [36|. This is not feasible, in 
practice, given the sizes L = 10 7 we study here. 

Instead, we used the circulant embedding method for 
fast generation of Gaussian field proposed in [2(| [2l| . 
This method allows to generate random increments which 
are approximately correlated according to Eq. ([2]) by gen- 
crating a periodic increment sequence of period L' , with 
L' > 2L. The correlation of this periodic sequence are 
encoded in a covariance matrix C;.;+ m = C(m) of size 
L' x L' built using the original covariance C and defined 
as: 

C(m) = C(m) for m = 0, . . . , L'/2 - 1 (3) 
C(m) = C(L' - m) for m = L'/2, L' — 1 . 

Toeplitz matrices displaying this periodicity are called 
circulant matrices. For the actual analysis of the nu- 
merical simulations we consider only the first L steps 
Ax\, A.X2, ■ ■ ■ , Axl- If L is large, the correlation between 
the first and the last increment is small and the periodic- 
ity has no large influence. The advantage of this approach 
is that the periodicity of the matrix C allows the appli- 
cation of Fast Fourier Transformation (FFT) to generate 
the fflms [12. 

The FFT is performed in ~ L'log(L') operations, 
technically we use the GNU Scientific Library (GSL) 
[H[. Let be c k the FFT of C(m) from Eq. i.e., 

Cfe = J2m=o C(m)e~ 27 "'~ m . Since C(m) is symmetric 
and positive, the coefficients cu are real positive num- 
bers. The generation of the correlated random numbers 
works as follows. 

(i) the starting point are L' independent and identically 
distributed (iid) Gaussian numbers (with zero mean and 
variance one). 

6 = (£oj£i) • 

(ii) we define 

h = y/VZk£k (fc = 0,...,£'-l) (4) 

which are real numbers as well and where the factor \JU 
takes into account the correct normalization. 

(iii) The vector of the increments is obtained after back 
transforming the vectors of Sk- 

5i = ^Y^he 2 ^\ (5) 

k=0 
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and taking real and imaginary part: 

Ax t =Re{^} + Im{^} (6) 

with I = 1,2,..., L'. 

It is easy to check that these three steps lead to the 
desired correlation. Using = and — Sk,k' 

one arrives at (5xf^ 5x^) = C(±i ± j) where the first 
and second signs are + for the case of having Sxi and 
Sxj on the left, respectively, and — for the conjugate 
complex, respectively Thus, using Re(z) = (z + z*)/2, 
Im(z) = (z — z*)/2i, and C(m) = C(— m) one obtains 
finally, as desired 

(Axi +m Axi) = 

-((Si +m + 5* l+m - iSi +m + iSl +m )(Si + 6* - iSi + iS*)) 

= C(m) . 

A numerical test of the method is shown in Fig. [TJ which 
proves that indeed the generated ranomd numbers follow 
Eq. ©. 




m 



FIG. 1: Correlation between increments. Shown are the 
wanted function C(m) and the numerical data for L — 10 4 . 
The main plot is for H = 1/4, which the inset displays 
H = 2/3. 

For a direct simulation of fBm processes, one gener- 
ates a vector of V real random numbers, constructs the 
vector of complex numbers, £, uses the transformation 
Eqs. (HI,® to obtain the correlated increments Axi and 
finally 

L 

x = x(L) = Ax i ( 7 ) 
i=i 



Nevertheless, since we use an absorbing boundary at 
x = 0, most of the time at least one of the intermediate 
steps will visit the negative half axis, i.e. Y2i=i ^- x i < 
(for some £ < L), and the obtained value x(L) does not 
contribute to the distribution P + (y). The probability 
of being not absorbed, i.e., the persistence (or survival 
probability), behaves like 

S(x = 0,L) ~ L-° (8) 

with 9 being t he p ersistence exponent known to be 9 = 
1 - H 0, [H HI- Hence, for the case L = 10 7 and 
H = 1/4, which we study here, we obtain Po(L) ~ 10~ 5 . 
This means, a direct simulation is not feasible. The direct 
simulation approach has been used in the past [l9[ only 
for values of H > 1/2, which is simpler than H < 1/2. 
Nevertheless, due to the mentioned limitations, even for 
the simpler case only sizes of L = 20000 could be studied, 
in contrast to L = 10 7 , which is studied here for values 
H > 1/2 as well as for H < 1/2. Note that an alternative 
is to directly simulate a physical process which exibits 
the nature of fractional Brownian motion, like a suitable 
polymer model. Nevertheless, this includes many more 
physical details as a raw fBm, hence only even smaller 
system can be studied, like in a recent study [2a |. where 
polymers of length N = 257 monomers could be treated. 

To circumvent this problem, we performed a Markov- 
chain Monte Carlo simulation with the configuration 
space being the set of all feasible random vectors £. Fea- 
sible means that the resulting fBm process (after FFT 
to generate the correlation of the increments Ax) is not 
absorbed. The simulation must be initialized with an al- 
lowed configuration, namely we start from a random vec- 
tor and a corresponding correlated increment Aa;*- ^ 
such that the resulting process is not absorbed. In prac- 
tice, for H < 1/2, we facilitate the generation of a feasible 
initial configuration by sampling from a shifted Gaussian 
1) (£ > 0) and repeat the search for an initial con- 
figurations until one feasible increment vector is found. 
This initial configuration is clearly biased, but does not 
have influence on the final result since only after some 
sufficient equilibration time we start to sample the ob- 
servables. 

Each Monte Carlo step £^ — > consists of chang- 

ing a fraction p of randomly chosen entries of the config- 
uration fW, the new entries being iid G(0,1), resulting in 
a trial configuration £ tIlal . Then, again after using FFT 
to introduce the correlation, we obtain Ax tria . if the re- 
sulting fBm {J2f = i Axf ia1 } (L = 1, . ..,L) is absorbed, 
the trial configuration is rejected, i.e., £( t+1 ) = £w. If 
the resulting fBm is allowed, the trial configuration is 
feasible, hence it is accepted, i.e., £( t+1 ) = £ trial . This 
approach satisfies detailed balance, hence converges to 
the correct distribution: The distribution of the configu- 
rations is given by a product of Gaussians over the space 
of feasible configurations, i.e., 

P(0 = n(^=exp(-4 2 /2))/« 
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where the indicator function 1^ is 1 if £ is feasible, i.e., 
the resulting fBm is not absorbed, and else. Hence, if 
a certain fraction p of the entries of £ is replaced to yield 
the resulting change of weight is given by 



1 



exp(-ef /2) h 



where the product runs over the changed entries. This 
change of weight is symmetric to the exchange £ -H- 
hence detailed balance is fulfilled: P(£,) w (£ ~^ £') — 

The Markov chain in the configuration space 
is reflected by the sequence of the endpoints 
x^-°\L) , x^\L) , . . . of our Monte Carlo simulation. 
Here we studied the statistics of the rescaled variable 
y = x(L)/ \<jL h ) Since we are interested in the behavior 
of P+(y) near y — 0, we also used a bias b(y) = y~ a 
(a > 0), by imposing an additional Metropolis criterion 
[27H29| and accepting a feasible configuration with 
the probability p accopt = min{l, 6(t/ w )/6(j/W)}. This 
drives the simulation into the range of interest. We 
adjusted the fraction p of changed entries such that the 
total acceptance probability of an MC step is near 0.5. 
Hence, for each value H of the Hurst exponent and each 
length L, we had to find a suitable value p = p(H, L) 
empirically. 
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t 

FIG. 2: Sample trajectory of a Monte Carlo simulation: End- 
point y (t) = x{L) {t) /L H of a non-absorbed fBm (H = 1/2) of 
length L — 10 7 as a function of the MC time t, for the initial 
phase of the MC simulation up to t = 10 4 . 
used to increase the statistics near y — 0. 



A bias 



In Fig. [2] a sample trajectory in the space of (rescaled) 
endpoints is shown for the case H — 1/2 and L = 10 7 . 
Via a bias b(y) ~ y~ 12 the simulation is concentrated 
near y sa 0. 
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H = 1/4 


H = 2/5 


H = 1/2 


H = 2/3 


H = 3/4 


10 3 


0.010 


0.060 


0.030 


0.40 


0.50 


10 4 


0.020 


0.020 


0.020 


0.40 


0.50 


10 s 


0.020 


0.030 


0.020 


0.20 




10 6 


0.010 


0.020 


0.020 


0.10 




10 7 


0.015 


0.020 


0.015 


0.05 


0.05 



TABLE I: Value of the Monte Carlo parameter p for different 
lengths of the fBms and different values of the Hurst exponent 
H. 



Concerning equilibration of our Monte Carlo Simula- 
tion, we found that typically, for the longest fBm pro- 
cesses, after 1000 sweeps we do not find any sign of the 
initial configuration. After disregarding this initial bunch 
of Monte Carlo sweeps, we measured histograms [3(| of 
the rescaled endpoints of the processes. In case a bias is 
applied, the histograms have to be multiplied by = y a 
and normalized to get the final distributions P+{y). 

Note that the approach use here is rather general: Dur- 
ing the Monte Carlo simulation a vector of variables is 
changed. This vector is evaluated and within a Metropo- 
lis criterion is is determined whether the changed vector 
is accepted. This is like in any Markov chain Monte Carlo 
simulation, e.g., a single-spin flip simulation of the Ising 
model. The only difference is that the step transform- 
ing the configuration vector into a Metropolis cirtierion 
is rather involved here since it includes creating a corre- 
lation between the vector entries, turning them into ran- 
dom walks, checking for absorption and including a bias 
keeping the random walks close to the origing. For the 
Ising system the same step would be just the calculation 
of an energy. Nevertheless, this approach allows to treat 
non-equilibrium non-stationary processes, like fractional 
Brownian motion, within the same framework as conven- 
tional equilibrium statistical mechanics systems. Hence, 
the approach should be applicable to a wide range of 
problems. 



III. RESULTS 

We have performed simulations to generate fBms for 
values of the Hurst exponent H = 1/4, H = 2/5, 
H = 1/2, H = 2/3 and H = 3/4 of lengths L = 10 3 , 
10 4 , 10 5 , 10 6 and 10 7 , respectively (for H = 3/4 we did 
not consider L = 10 5 and L = 10 6 since this was not nec- 
essary). For the rescaling, we used a = 2 (H = 1/4), 
a = 1.5 (H = 2/5), a = 1.2 (H = 1/2), a = 0.5 
(H = 2/3), and a = 0.4 (H = 3/4), respectively. For 
each case, we determined the parameter p, such that the 
acceptance probability is (very roughly) about 0.5. The 
values we used are shown in tableQ] For each case the MC 
simulation was performed for long runs, up to t — 3 x 10 8 
for the longest walks of length L — 10 7 . 

Note that for H = 1/4, H = 2/5, and H = 1/2, we 
have restricted the simulations to fBm processes with y > 
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0.0001 to prevent the simulation being caught near y = 
due to a very small acceptance ration via the rescaling 
factor in that region. 
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FIG. 3: Distribution P+(y) of rescaled endpoints y = x/L H 
for non-absorbed fBms (Hurst exponent H — 1/2) 

First, to verify our method, we studied the case of 
standard random walks, H = 1/2. For this case it is 
possible to know the first corrections to the continuum 
limit behavior (see Appendix): 



P + (y,L) = f (y)--j=f 1 (y) 



(9) 



where L is the number of increments, the constant c de- 
pends on the increments distribution of the random walk. 
For Gaussian numbers (zero mean, unit variance) we have 
c = C(1/2)/v / 2tt 0.582597. . . and the scaling func- 
tions are: 



fo(y) = ye 



Mv) 



(10) 
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V/2 



The rescaled distributions P+{y) of the endpoints are 
shown, together with the predictions of Eq. (|9j) valid for 
large L, in Fig. [31 One is able to see strong finite-size 
effects for small values y — > 0, where a plateau is visible. 
For increasing length L, the plateau decreases as c/vL 
and the data approaches better and better the continuum 
limit scaling function /o. 

We conclude that for a generic fBm for which first cor- 
rections to the continuum limit behavior is not known, 
the plateau should also vanish when the size of the system 
is large and for very long walks (L — 10 7 ) the continuum 
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FIG. 4: (color online) Distribution P+{y) of rescaled end- 
points y — x/L H for non-absorbed fBms (Hurst exponent 
H = 2/3) 



limit behavior is displayed over several order of magni- 
tudes. 

Based on scaling arguments it has been conjectured Q 
that, in the continuum limit P+(y) vanishes as y^ with 
= (1 — H)/H = 9/H for y — >• , 6 being the persis- 
tence defined via Eq. ([8]) . This conjecture was confirmed 
by an epsilon expansion around the Brownian solution 
obtained thanks to a field theory calculation [l9j. The 
numerical check of the conjecture for values of H far from 
1/2 remains very challenging. We first consider the dis- 
crete random walk with H = 2/3. In this case, since the 
persistence is decreasing not very fast (9 — 1/3), numeri- 
cal results were obtained for moderate lengths L — 2 x 10 4 
by direct simulations fl9l ]. which were compatible with 
the analytics. Here, we were able to study this case again. 
Our results, up to a length of L = 10 7 , confirm the ana- 
lytics with much better accuracy, as the small-endpoint 
behavior follows the expected power law with exponent 
0=1/2 very well, see Fig. 0] Also for the case H = 3/4, 
the behavior close to the origin matches the expected 
P{y) ~ y^ behavior with = 1/3 very well, see Fig. [5] 

For the subdiffusive case the convergence with the 
length L of the path is slower. We first consider the case 
H = 2/5, see Fig. [S] For the longest walk, the behav- 
ior close to the origin follows a power law P(y) ~ y^" 11 , 
but the exponent cf> c s — 1.44(1) (obtained from fitting a 
power law in the region y E [0.03, 0.2]) is slightly smaller 
than the predicted value (j> = (1 — H)/H = 3/2. A bet- 
ter estimation is obtained by performing a finite-length 
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FIG. 5: (color online) Distribution P+(y) of rescaled end- 
points y — x/L H for non-absorbed fBms (Hurst exponent 
H = 3/4). 
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FIG. 6: (color online) Distribution P+(y) of rescaled end- 
points y — x/L H for non-absorbed fBms (Hurst exponent 
H = 2/5) 



extraplation of the form 

4> cS {L) = <b + cL- b (11) 

for the effective exponent as a function of the length L, 
see Fig. [7] When fitting Eq. (fTT|) to the data, we obtained 
4> = 1.50(4) (and b = 0.23(4)), in prefect agreement with 
the prediction. 

Finally, we turn to the most difficult case whith H = 
1/4 where we expect that P+(y) vanishes as y 3 as y — > 0. 
Direct simulations on this process (for restricted sizes) 
are not conclusive and a scaling behavior ~ y 2 was found 
to be consistent with the data [2(| [3l[ . Using our Markov 
chain approach we can see the finite size effects remain 
important even for long processes (see Fig. [SJ. Again, 
like in the case H = 2/5, we observed power-law behv- 
ior close to the origin with an effective exponent cj) c s(L). 
From L — 10 5 on, the effective exponent is above 2 and 
is growing with walk length L. Hence, previous claim 
of a quadratic behavior can be discarded clearly. Here 
the extrapolation according to Eq. (fTTj) yields <f> = 3.3(3) 
with a very slow convergence b — 0.08(2). Hence, even if 
one was able to study extremely long walks like L = 10 11 , 
one would observe </> e ff ~ 2.88. Thus, the observed ex- 
trapolated exponent is also compatible within error bars 
with the predicted value = 3, but with lower accuracy 
due to stronger finite-length correction. 



IV. SUMMARY 



We have introduced a Markov-chain Monte Carlo ap- 
proach to study numerically fractional Brownian motion 
in the presence of an absorbing boundary via generating 
finite-step random walks with correlated disorder. Our 
approach allowed us to study long walks up to L = 10' 
steps. For the test case H = 1/2 the result for the dis- 
tribution P+{y) of the rescaled endpoints y of the walks 
agrees in the limit L — > oo with the exacty analytic result. 
We also derived analytical expression for the fite-length 
corrections, which turn also to be compatible with the 
numerical results, better with increasing step number L. 

In the main part of our work, we studied fractional 
Brownian motion where we find for y — > power-law 
behaviors P+{y) ~ yr . For the superdiffusive cases, 
H = 2/3 and H = 3/4, we observed for long walks 
L = 10 that the measured exponents match the analyt- 
ical prediction tj> = (1 — H)/H with very good accuracy. 

For the subdiffusives cases H — 2/5 and H = 1/4 we 
found strong finite-length effects which can be described 
via an effective exponent (f> e g(L). Hence we could not 
observe the limiting exponent directly, but in both cases 
we found via a power-law extrapolation a convergence to 
the predicted values <f> = (1 — H)/H. 
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FIG. 7: Effective exponent e g as a function of the walk 
length L for // = 3/4. The line shows a fit to the function 
<f>es(L) = 4> + aL~ b . Inset: for the case H = 1/4. 
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FIG. 8: (color online) Distribution P+{y) of rescaled end- 
points y — x/L H for non-absorbed fBms (Hurst exponent 
H = 1/4) 
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Appendix A: Derivation of Eq. 8 and Eq. 9 



density that the particle arrives at x at step L while 
staying above at all intermediate steps. An exact ex- 
pression for pl(x), or rather for its generating function, 
is known explicitly for arbitrary jump density f(i]) and 
is given by [321 ] 



/ dx e~ x x V p L (x) s L = (j>(s, A) 



with 



^( S ,A)=ex P --/ fc2 + A2 dk 



(A2) 



(A3) 



where f(k) = e lkv f(rj) drj is the Fourier transform 
of the noise density. Our goal is to extract the leading 
(and subleading) scaling behavior of pl(x) for large L 
from Eq. (f&2|) . 

To make progress, it is useful to consider an alternative 
expression for <j>(s,\) derived in Ref. [33|, valid for all 
/(f?)'s with a finite variance tr 2 = rf f(rj) drj, 



We consider a random walk starting at the origin. Its 
position at discrete time steps evolves via 



•En — 1 4" 



(Al) 



starting from xq = 0. The random variables rj^s are in- 
dependent and identically distributed noises, each drawn 
from a symmetric and continuous probability density 
function (pdf) f(r)). Let pl{x) denote the probability 



(s,A) 



x exp 




We next consider the scaling limit when x — !• oo, L —> oo, 
with the ratio y = x/^/L fixed. In the Laplace place, this 



8 



corresponds to taking the limit A — > 0, s — > 1, keeping 
the ratio X/y/1 — s fixed. Taking this scaling limit in Eq. 
(|A4|) . one gets 



1 



•A 



vT^ + ctA/x/2 



(A5) 



where c is a constant with the following expression [33j, 

HI 



1 

c = — 

7T 



i - /(fc) 

<r 2 fc 2 /2 



(A6) 



Substituting the scaling-limit expression of cf>(s, A) from 
Eq. (IA5P on the right hand side of Eq. (| A2|) and inverting 
the Laplace transform with respect to A gives, 



L=0 



V2 



i + — cVT 



,-^2(l-s)x/o 



(A7) 

valid in the scaling limit s — > 1, x — > oo but keeping the 
product yl — s a; fixed. Next, one can invert this gener- 
ating function with respect to s, using Cauchy's inversion 
formula. Skipping details, we find that the two leading 
terms, in the scaling limit where x —¥ oo, L —¥ oo, but 
keeping y — x/yL fixed, are given by 



Pl{x) 



ye 



-y /2(7 _ 



7T 



2 /o 2 



(A8) 



The conditional probability Pl (x) (probability density 
to reach the position x given that it has survived up to 
L steps) is defined as 



Pl(x) 



Vl{x) 
f™PL(x)dx 



(A9) 



Substituting the scaling behavior for pl(x) from Eq. 
(|A8p in the above definition, we find that Pl{x) has the 
following scaling behavior 



P L (x)^^=P + (y,L) 



(A10) 



with y — xj \J~L and 

P+(y,L) = f (y) 
where 



VI 



h(y)+0(l/L) 



My) = 
My) 



= JL P -v 2 /2° 2 



a- 

e~ 



2 / 2 2 2 / n 2 

v y e -v /2o- 

7TCT 



(All) 

(A12) 
(A13) 



and the constant c is given by Eq. (|A6|) . For the special 
case of the Gaussian jump density, f(?]) = e~ v I 2 jsf^K 
(with er 2 = 1), one can evaluate the constant c in Eq. 
HI explicitly H3 



C 1/2 

c = sv /_ ; = -0.582597. 



In this case, in particular, putting y = we get 

0.582597... 



(A14) 



P + (0,.),-- = 



(A15) 



which is consistent with our simulations. 
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